Transport in anisotropic model systems analyzed by a correlated projection 

superoperator technique 
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By using a correlated projection operator, the time-convolutionless (TCL) method to derive a 
quantum master equation can be utihzed to investigate the transport behavior of quantum systems 
as well. Here, we analyze a three-dimensional anisotropic quantum model system according to 
this technique. The system consists of Heisenberg coupled two-level systems in one direction and 
weak random interactions in all other ones. Depending on the partition chosen, we obtain ballistic 
behavior along the chains and normal transport in the perpendicular direction. These results are 
perfectly confirmed by the numerical solution of the full time-dependent Schrodinger equation. 
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I. INTRODUCTION 

The transport of different extensive quantities like en- 
ergy, heat, entropy, mass, charge, magnetization, etc., 
through and within solid state systems is an inten- 
sively studied topic of nonequilibrium statistical dynam- 
ics. Nevertheless, there are numerous open questions con- 
cerning the type of transport especially in small systems 
far from the thermodynamic limit and in particular in 
quantum mechanics. At the heart of many investigations 
is the classification into two main categories: normal or 
diffusive transfer of the extensive quantity and ballistic 
transport featuring a divergence of the conductivity. 

Diffusive transport occurs whenever the system is gov- 
erned by a diffusion equation. In particular, this means 
that excitations decay exponentially fast and the spa- 
tial variance of an initial excitation grows linear in time. 
Ballistic transport, however, is rather described by the 
equations of motion of a free particle. For the spatial 
variance of an excitation this implies a quadratic growth 
in time. 

In the present paper, we will concentrate on the trans- 
port of energy and heat in quantum systems. There are 
several different approaches discussed in the literature to 
investigate the transport of those quantities in quantum 
mechanics. One very famous ansatz is the investigation 
of heat transport in terms of the Green-Kubo formula 
[1-6]. A main advantage of this approach is certainly 
its computability after having diagonalized the system's 
Hamiltonian. Derived on the basis of linear response the- 
ory the Kubo formula has originally been formulated for 
electrical transport [7, 8], where an external potential 
can be written as an addend to the Hamiltonian of the 
system. Basically one finds a current-current autocorre- 
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lation, which has ad hoc been transferred to heat trans- 
port simply by replacing the electrical current by a heat 
current [9]. However, the justification of this replace- 
ment remains unclear since there is no way of express- 
ing a temperature gradient in terms of an addend to the 
Hamiltonian of the system as before [10]. 

Other approaches to heat conductivity in quantum sys- 
tems are based on diagonalization of the Schrodinger 
equation [11], analyzing the level statistics of the Hamil- 
tonian [12, 13] or by an explicit coupling to some envi- 
ronments of different temperature [14, 15]. In the latter 
case, environments are described by a quantum master 
equation [16] in Liouville space. Here the temperature 
differences can, indeed, be described by a perturbation 
operator so that one may treat a thermal perturbation 
in this extended state space similar as an electrical one 
in the Hilbert space [17]. 

The Hilbert space Average Method [18] allows for a di- 
rect investigation of the heat transport in quantum sys- 
tems from Schrodinger dynamics. By deriving a reduced 
dynamical equation for a class of design quantum sys- 
tems, normal heat transport as well as Fourier's Law has 
been confirmed [19, 20]. Recently, it has been shown that 
for diffusive systems the Hilbert space average method is 
equivalent to a projection operator technique with an ex- 
tended projection superoperator [21, 22]. However, bal- 
listic behavior cannot be analyzed with the Hilbert space 
Average Method in a straightforward manner since it is 
not obvious how to obtain time-dependent rates. 

Using a correlated superprojection operator within the 
derivation of the time-convolutionless (TCL) quantum 
master equation leads to a reduced dynamical descrip- 
tion of the investigated system. The main advantage of 
the correlated TCL method refers to its perturbation the- 
oretical character. Thus it is a systematic expansion in 
some perturbational parameter. 

To use this alternative method for an investigation of 
the transport behavior of a quantum system, it is nec- 
essary to partition the microscopic system described by 
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the Hamiltonian H into mesoscopic subunits. While the 
complete dynamics is governed by the Schrodinger equa- 
tion of the full system according to its density operator 



[H,p\ = C{t)p, 



(1) 



wo aim at deriving a closed reduced dynamical equation 
for the subunits chosen. Formally, this partitioning is 
done by introducing a projection superoperator V that 
projects onto the relevant part of the full density matrix 
p [16], here the spatial energy distribution within the 
system. However, by a straightforward application of 
the projection superoperator on the above equation, the 
dynamics of the reduced system is no longer unitary, but 
described by 



Vp = VC{t)p. 



(2) 
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FIG. 1: Partition schemes for investigating tiie transport per- 
pendicular (a) or parallel (b) to the spin chains. Each spin 
is represented by a dot, solid lines indicate Heisenberg in- 
teractions along the chains, dashed lines represent random 
interactions. The diagonal couplings within each plane have 
been left out for clarity [except the lower left corner of (a)]. 



These effective equations of motion for the relevant 
part Vp can either be written as an integro-differential 
equation (Nakajima-Zwanzig equation [23, 24]) or as a 
time-convolutionless (TCL) master equation [16], which 
is an ordinary linear differential equation of first order. 
Both methods allow for a systematic perturbative expan- 
sion. In the TCL expansion series the first-order term 
typically vanishes and thus the leading order is given by 
(cf. [21,' 22]) 

t 

p^ = - j dtiV£{t)£{ti)Vp. (3) 



However, in order to obtain a converging perturbation 
series expansion V should not be chosen arbitrarily: A 
"wrong" projection superoperator may lead to a break- 
down of the expansion [21]. 

II. DESCRIPTION OF THE MODEL 

In the present paper we consider a three-dimensional 
(3D) model composed of two-level systems. The coupling 
between the atoms is anisotropic, i.e., in one direction 
dominated by a Heisenberg interaction whereas the cou- 
pling in all other directions is random. The choice of 
random couplings ensures that the interaction is unbi- 
ased as it does not have any special symmetry. Two- 
level atoms or spin-1/2 systems [25] allow to study a 
large variety of quantum effects from quantum infor- 
mation processing to solid state theory, described by a 
rather simple interaction, making them interesting both 
from an experimental and theoretical point of view. Of 
particular interest are the transport properties of sys- 
tems containing ID and 2D spin structures, e.g., the in- 
vestigations of heat transport in cuprates, in which a 
dramatic heat transport anisotropy has been reported 
[26, 27]. While the anisotropy is mainly attributed to 
anisotropic phonon scattering processes, we investigate 



transport anisotropics emerging from an anisotropic (but 
coherent) interaction. 

The model we are going to investigate is a three- 
dimensional model of two-level systems depicted in Fig. 1. 
In terms of Pauli operators the local Hamiltonian of the 
network is given by 

Hloc = ^T.-^\ (4) 

i 

with the local energy splitting AE defining the basic en- 
ergy unit within our model. 

In X direction, the two-level systems are coupled via a 
Heisenberg interaction 

^^ = ^6-«®<t('+^\ (5) 

i 

with the Pauli spin vectors (T^*^ = {c7x \(yy\(7i^^) at site 
i. 

In the y and z directions, we use a random interac- 
tion matrix to couple both adjacent sites and next 
neighbor sites lying diagonally opposite [sec lower left 
corner of Fig. 1(a)]. The nonzero matrix elements are 
taken from a Gaussian ensemble with zero mean and a 
variance s^. While each matrix element is taken from 
the same ensemble, the geometry of the system requires 
that we do not have translational invariance within the 
random interaction. 

To investigate the transport in a; or 2 direction, re- 
spectively (cf. Fig. 1), we perform a partition of the 
model into N subunits. A layer of n two-level systems 
is grouped together into a new local subsystem, coupled 
to adjacent layers by the connections between pairs of 
two level systems. Because of the anisotropy within the 
model we can study the transport perpendicular to the 
Heisenberg chains in the z direction [Fig. 1(a)] and along 
the chains in x direction [Fig. 1(b)]. 
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FIG. 2: N subunits with ground state and first excitation 
band of width 6s containing n energy levels each. Black dots 
specify the initial states used. 



The coupling strength of an arbitrary interaction ma- 
trix V is defined as 



(6) 



with d being the dimension of the matrix (see [18]). For 
the random interaction we choose the variance in such 
a way that 77 = 1 for all interaction matrices coupling 
adjacent subunits. 

The complete Hamiltonian of the full model system is 
thus described by 



H = Hioc + XhHh + XrHh. 



(7) 



Because of the normalization of the interaction matrices 
the numbers Xh and define the coupling strength be- 
tween different sites. The coupling strengths Xh for the 
Heisenberg interaction and for the random interac- 
tion are chosen so that Xr <C Xh <C AE, which is known 
as the weak coupling limit. 

Regardless of the partition scheme chosen (in the x ov z 
direction) each subunit can be seen as a molecule consist- 
ing of several energy bands. However, the solution for the 
complete system is computationally unfeasible for more 
than a few sites. If we restrict ourselves to initial states 
where only one site is excited (or superpositions thereof) 
the Heisenberg interaction does not allow to leave this 
subspace of the total Hilbert space. By choosing also 
the random interaction to conserve this subspace we re- 
strict all further investigations to the single excitation 
subspace. Figure 2 gives a graphical representation of 
our system, with 5e being the width of the first energy 
band (all higher excitation bands are neglected here). 

Interpreting our model system in terms of a magnetic 
system, i.e., the two- level atoms representing coupled 
spins in a magnetic field for example, the considered en- 
ergy transport is equivalent to spin transport in a gapless 
system (i.e. AE = 0). 



III. TRANSPORT IN THE z DIRECTION 

A. Partitioning scheme 

Let us consider the transport perpendicular to the 
Heisenberg chains (in the z direction) first. Then, the 



partitioning into subunits yields the following mesoscopic 
Hamiltonian consisting of a local and an interaction part 



H = Hl+Hj 

N N-1 

= ^^i(/x)+^^,(M,/x + l). 



(8) 



1^=1 



Here Hl{^) of subunit /j. consists of the constant local 
energy splitting, the Heisenberg interaction, and the in- 
ternal random couplings of each subunit [cf. gray planes 
in Fig 1(a)]. Since X^ <C Xh the effect of the internal 
random couplings on the spectrum of Hl{^x) may be ne- 
glected. Therefore the bandwidth 6e is determined by 
the Heisenberg interaction given by 



5e = 8A 



H- 



(9) 



The last term in Eq. (8), Hi{jjL,ii + 1), denotes the in- 
teraction between the subunits which is purely random 
here, i.e., contains parts of the random interaction Hamil- 
tonian Hr only. 



B. Derivation of the TCL master equation 

The correlated projection superoperator V introduced 
in Sec. I is of the type as suggested by Breuer [22] and 
reads 



Pp = ^Tr{n,p}ln,^5]p,ln„ 



(10) 



with being the standard projection operators 



and |n^) the eigenstate of Hl{i^) in the one-particle ex- 
citation subspace, i.e. the states in the band of subimit 
/Lt (cf. Fig. 2). Consequently, the number is just the 
excitation probability of subunit /U. This choice of V thus 
implements the partitioning scheme required for studying 
transport behavior. 

Switching to the interaction picture, plugging both the 
Hamiltonian (8) and the projection (10) into Eq. (3) we 
get 



X2 



dtiTr{n^[^K(t), \HR{ty)X\\}P. 

(12) 

for the second order TCL expansion. The time dependen- 
cies of the coupling operators refer to the transformation 
into the interaction picture and are defined as 



(13) 
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By exploiting that 11^ projects onto eigenstates of 
Hl (a*) we can evaluate the trace by using the block struc- 
ture of the interaction Hi{ij., /x-I- 1) between adjacent sub- 
units (see [19, 20]), resulting in 



ii = -7^(2P^-P^+i-P^_i) 



with the decay rate 



2A 



1, = ^T.\^K\Hr\1,+i)\- 



sm{LUkit) 



k,l 



(14) 



(15) 



The frequency w^; refers to the transition between the 
eigenstates k, I. Equation (14) is basically a rate equation 
for the probabilities to find an excitation in subunit ji. 



C. Decay rate 

Since the interaction between two adjacent subunits 
is a random matrix with the above described proper- 
ties, all matrix elements are approximately of the same 
size. That means that the rate does not depend on 
the subunit ji (7^ = 7). Furthermore, we can assume 
|(fc^|7?fl|;^+i)|2«l, finding 



(16) 



In the following the double sum is treated analogous to 
the derivation of Fermi's Golden Rule. 

Since the sine cardinal (sine) of Eq. (16) is a represen- 
tation of the Dirac (5-distribution 



7i"5t(wfe0 = lim 



sin(wfcit) 



we may approximate the rate for not too small t by 



7; 



nh 



(17) 



(18) 



Replacing the double sum over integrals in the energy 
space we arrive at 



1' 



9^{E)dE 



(19) 



with the state density g{E), i.e., the integral over the 
square of the density of states. Since we have neglected 
the internal random interaction completely, the state 
density of the first excitation subspace is just given by 
the state density of a Heisenberg spin chain 



9{E) 



2n 



1 



(20) 



Unfortunately, this function is not square integrablc due 
to singularities at the boundaries of the spectrum. How- 
ever, due to symmetry we have 



/ ' g^{E) dE = 2 r' g\E) dE. (21) 
Jo Jo 

In order to avoid the singularity at £J = we renormal- 
ize the number of states in the band. We introduce the 
regularized integral 

/.(5e/2 

FA{n) = 2 a''g''{E)dE 

J A 

with a being the factor that renormalizes the number 
of states. We assume that for a band consisting of only 
a few levels n (but still enough to define a density of 
states), the density of states is approximately constant. 
For a constant density of states g{E) we simply have 



r5e/2 

Jo 



f{E)dE='j-, 



(23) 



therefore our renormalization prescription is given by 



FK{n) = 



5£' 



(24) 



Using this result to solve Eq. (22) for a at constant h 
yields 



(25) 



V'21n(fc/A- 1) 



This allows us to calculate the physical limit of the renor- 
malization procedure, i.e.. 



limF4n) = - 



(26) 



which is the same value as for a constant density of states. 
This finally leads to the relaxation rate 
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(27) 



The approximation introduced by Fermi's Golden Rule 
is only valid in the linear regime (see [18]), i.e.. 



2„\2 
R 



fc2 



< 1. 



(28) 



D. Solution of the TCL master equation 



Figure 3 shows both the numerical results for the so- 
lution of the full Schrodinger equation and the solution 
of the rate equation (14), according to the above derived 
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FIG. 3: Perpendicular transport: probability to find the ex- 
citation in subunit ^ = 1,2,3. Comparison of the numerical 
solution of the Schrodinger equation (crosses) and second- 
order TCL (lines). (iV = 3, n = 600, Xr = 5 ■ 10"* AS, 
Ah = 6.25 ■ 10"'^ AE) 



the interaction transforms into 

= e'"'^' HRe-'"^\ (32) 

where Eq. (31) has been used. Note that this is not the 
standard interaction picture as used above, but a special 
one allowing us to treat the transport in the x direction in 
a similar manner as in the z direction. According to this 
transformation the derivation of the second order TCL 
master equation in Sec. IIIB, especially Eqs. (12), (14) 
and (15), remain imchanged. 

B. Decay rate 



approximation for the rate 7 [cf. Eq. (27)]. Both are in 
reasonably good agreement. 

Equation (14) is a discrete version of the diffusion equa- 
tion, which does not change when regarding the thermo- 
dynamic limit {n, N ^ 00, nXj^ = const). For a 5-shaped 
excitation at t = its solution is a Gaussian function, the 
variance of which grows linear in time. Therefore, it is 
evident that the heat transport is normal perpendicular 
to the chains. 



IV. TRANSPORT IN THE x DIRECTION 

A. Partitioning scheme 

In the following let us concentrate on the transport 
in the x direction, i.e., parallel to the chains. Thus, we 
have a slightly different partition of the total Hamilto- 
nian. Besides the local energy splitting, the local part 
Hl of the mesoscopic Hamiltonian (8) contains random 
interactions only: 

N 

Hl = Y. [Hloc{^i) + Afl7?fl(/i)] . (29) 

In contrast, the interaction between the subunits consists 
of a Heisenberg and a random part, 

N-l 

Hi=J2 [AHi?H(M,M+l) + Afli/fl(M,M+l)]- (30) 

In the one-particle excitation subspace the commuta- 
tor relations 

[Hh,Hl] = [Hh,Hu]=0 (31) 

are satisfied. If the dynamics induced by the local part 
Hl and the Heisenberg Hh is absorbed in the transfor- 
mation into the interaction picture, the random part of 



However, the computation of the rate (16) is different 
here. For calculating the local band structure we con- 
sider just a random matrix of dimension n, drawn from a 
Gaussian unitary ensemble. From random matrix theory 
[28] it is known that the density of levels C{^) for such 
a random Hermitian matrix consisting of elements with 
zero mean and unit variance for both real and imaginary 
parts is given by 



C{x) = -v/2n-a;2. (33) 
Mapping this to a density of energy levels leads to 

We rescale the variance to the interaction strength A^, 
which gives for the bandwidth 

6e = 4V^Afl. (35) 

In order to check whether our local Hamiltonian Hl 
can be approximated by such a random matrix, we com- 
pare the eigenvalues E{x) of both matrices. Using 



dx g[E{x)] ' ' 

and separating variables yields 

^J^jl-E^dE^dx, (37) 
noe V 4 

with the state density (34). This expression cannot be 
solved analytically for E, so we compare the numerical 
solution for discrete values of x with the eigenvalues of 
Hl- As Fig. 4 shows, Hl may indeed be approximated 
by a random matrix drawn from a Gaussian unitary en- 
semble. However, by plugging Eq. (35) into Eq. (28) one 
gets a constant value of 7r^/4 which is definitely not small 
compared to one. Thus, the requirement for the linear 
regime is violated and the derivation of the rate according 
to Fermi's Golden Rule can no longer be applied. 
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FIG. 4: Comparison of the eigenvalues of Hl and a random 
matrix drawn from a Gaussian unitary ensemble, (n — 600, 
Afl = 5 ■ 10-^AS) 
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FIG. 5: Parallel transport: probability to find the excitation 
in the first subunit (/i = 1). Comparison of the numerical 
solution of the Schrodinger equation (crosses) and second- 
order TCL (solid line). (Same parameters as for Fig. 3) 



The approximation used in Sec. Ill B is analogous to 
Fermi's Golden Rule. All transitions in Eq. (16) from I to 
k are weighted by the respective value of the sine function 
which changes its shape for increasing times to approach 
a delta peak for t — > oo. In the situation described above 
the decay takes place within the linear regime, i.e., at 
an intermediate time scale. That means that all possi- 
ble transition frequencies are distributed below the peak. 
Thus the sum in Eq. (16) can be approximated by the 
area under the peak (see [18]). 

This is not the case here. The decay happens on a 
much shorter time scale, when the peak is extremely 
broad. Therefore almost any transition frequency be- 
longs to the maximum of the peak. Thus the sine in 
Eq. (16) should better be approximated by the maxi- 
mum value of the peak, instead of the area under the 
peak. The maximum value grows with time according 
to t. Thus the double sum over sine functions could be 
approximated by n^t. This means that we get the relax- 
ation rate 

7=-^t (38) 



C. Solution of the TCL master equation 

The solution of Eq. (14) with the diffusion coefficient 
(38) defines the occupation probabilities in the interac- 
tion picture -P^"*. Note that in the other direction the 
occupation probabilities of the interaction picture have 
been equivalent to the occupation probabilities in the 
Schrodinger picture. This is not the case for the present 
situation because of the special choice of the interaction 
picture. Remember that we have used not only the local 
Hamiltonian for the transformation into the interaction 
picture, but also a part of the inter-subsystem interaction 
[cf. Eq. (32)]. 

Since we are interested in the occupation probabilities 
in the Schrodinger picture we need to calculate the 
inverse transformation of the density operator 

Vp- = e"""*Vp''''e'^"\ (39) 



where the diagonal elements "Pff^^ are the occupation 
probabilities P^. The off-diagonal elements of 7^/5'"* can 
be computed by replacing the projector (11) with another 
one projecting out off-diagonal elements as well. The 
dynamics of the diagonal and the off-diagonal elements 
decouple so that diagonal initial states remain diagonal 
for all time. 

Thus using Eq. (39) for the inverse transformation we 
get the time-dependent solution of the probabilities in the 
Schrodinger picture. In Fig. 5 the numerical solution of 
the Schrodinger equation is compared with the TCL pre- 
diction. Again, there is a very good agreement between 
the exact solution and our second order approximation. 



D. Spatial variance 

To classify the transport behavior in the x direction 
a very large system has to be considered, so that the 
initial excitation does not reach the boundaries of the 
system during the relaxation time. Since the solution of 
the time-dependent Schrodinger equation becomes un- 
feasible the second-order TCL prediction has been used 
for subsequent numerical integration. The variance of an 
excitation initially at /i = /iq, 

N 

shown in Fig. 6 grows quadratic in time, i.e., the trans- 
port is ballistic. Here we have considered a system with 
N — 300 subunits and an initial excitation at /ip — 150 
solving the TCL master equation. This is also valid in 
the thermodynamic limit as j{t) does not change. Nu- 
merical investigations show that the transport behavior 
is largely independent of 7(t). Ballistic transport is ob- 
served as long as Xnt ^ 7(t) on all relevant time-scales. 
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FIG. 6: Variance of an excitation initially at subunit no — 
150. Second-order TCL prediction (crosses) and quadratic 
fit (solid line). (iV = 300, n = 600, \r = 5 ■ 10"", 
\h = 6.25 • 10"^) 

V. CONCLUSIONS 

In the present paper we have demonstrated how the 
abstract method of correlated projection superoperators 
for the TCL master equation [21, 22] can be used to an- 
alyze the transport behavior of a three-dimensional solid 
state model: a system of coupled two-level atoms with 
an anisotropic interaction. The analysis is based on the 
following preconditions: 

1. a partitioning scheme in position space to consider 
the transport in one direction of the model, thus 
introducing a projection superoperator 

2. the convergence of the TCL expansion in second 
order (a wrong projection superoperator leads to 
a diverging expansion, or large higher than second 
orders) 

3. an approximation scheme for computing the decay 
rate to avoid numerical integration 

According to those central points a reduced dynamical 
description of the complex quantum model is derived 
which can be analyzed, e.g., to classify the transport be- 
havior of the system. 

By a comparison of the TCL prediction with the exact 
numerical solution of the complete Schrodingcr equation 
of our model system we have shown that the results of the 



method are in very good accordance with the real dynam- 
ical behavior of the system. Having established a method 
which efficiently describes the dynamical properties of a 
complex quantum model the transport behavior can be 
classified by either an analytic analysis of the solution 
of the reduced dynamical equations or by a numerical 
investigation. Here, the simplicity of the reduced equa- 
tions in comparison to the exact system of differential 
equations allows us to investigate the dynamical proper- 
ties of a much larger system which is not accessible to a 
direct investigation. 

The analysis shows that the model features two very 
different types of transport behavior in the x and z di- 
rections, perpendicular or parallel to the chains, respec- 
tively. In the z direction we have found a standard sta- 
tistical decay behavior following a diffusion equation on 
the basis of the mesoscopic subunits. In this way dif- 
fusive behavior has been derived from first principles 
on a mesoscopic scale whereas the dynamics on the mi- 
croscopic scale (i.e., of a single spin) is obviously non- 
diffusive. This indicates that the transport behavior is 
not only a property of a system per se, but also depends 
on the way we are looking at it. In contrast the model 
shows ballistic behavior parallel to the chains which is 
demonstrated by the features of the reduced dynamical 
equations. Note that this behavior is similar to investi- 
gations of large anisotropics within the heat conductivity 
of cuprates [26, 27]. 

In conclusion this method of a correlated projection 
superoperators within TCL allows to investigate the dy- 
namical behavior of 3D model systems on a mesoscopic 
scale. It is useful both in the case of a statistical decay 
according to a diffusion equation and the ballistic case, 
where time dependent rates are important. 
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